%
% Estimate productivity based on Ackerberg-De Loecker (2010)


% Get data
clear
workdir = 'E:/Dropbox/Work/Input_quality';
addpath = strcat(workdir,'/matlab');
path(addpath,path);

M = csvread(strcat(workdir,'/Data/matlab_data.raw'));  % Import CSV data file
M = sortrows(M,[9 1 2]);                          % Sort by nace, firm and year

% Generate some lags. 
M(:,12:16) = NaN(length(M),5);
  for i=2:length(M)                              
    if (M(i-1,1)==M(i,1) & M(i-1,2)==M(i,2)-1) 
      M(i,12)= M(i-1,4);  % Labor (quality adjusted)
      M(i,13)= M(i-1,10);  % Labor (unadjusted)     
      M(i,14)= M(i-1,5);  % Intermediates
      M(i,15)= M(i-1,7);  % Capital lag 1
    end
    if (i>=3)
      if (M(i-2,1)==M(i,1) & M(i-2,2)==M(i,2)-2) M(i,16)= M(i-2,5);  % Intermediates lag 2
      end    
    end
  end

% Load worker data
N = csvread(strcat(workdir,'/Data/matlab_data_wrk.raw'));  % Import CSV data file
N = sortrows(N,[3 1 2]);                                   % Sort by nace, firm and year

% Generate more lags. 
N(:,17) = NaN(length(N),1);
  for i=2:length(N)                              
    if (N(i-1,1)==N(i,1) & N(i-1,2)==N(i,2)-1) 
      N(i,17)= N(i-1,16);  % Average wage, t-1
    end
  end

nace = M(:,9);
tabulate(nace)

% Number of firm-years
disp('Number of firm-year obs'); 
disp(length(M));

% Number of firms
orgid = M(:,1);
tmp = unique(orgid);
disp('Number of firms'); 
disp(length(tmp));

% Number of firm-years with positive investment
posinv = M(:,6)>0;
disp('Number of firm-year obs with positive investment'); 
disp(length(M(posinv)));

% Number of firms with positive investment
tmp = unique(orgid(posinv));
disp('Number of firms with positive investment'); 
disp(length(tmp));

naceidx = [15 17 18 20 21 22 24 25 26 27 28 29 31 33 34 35 36 37];
J = 250;     % Number of bootstrap iterations
disp('Data loaded');

%% --------------------------------------------
% ESTIMATION : UNADJUSTED LABOR
% --------------------------------------------

export_unadj = [];

for nace = naceidx 
  disp('NACE'); disp(nace);
  firmidx = M(:,9)==nace;
  Mnace = M(firmidx,:);
  Nnace = N(firmidx,:);
  
  [beta export] = estimate_unadj([Mnace Nnace]);

  % Boostrap standard errors
  beta_bs = [];
  for j=1:J
    disp('Bootstrap iteration'); disp(j);
    Mboot = bootstrp2([Mnace Nnace]);                            
    beta_bs(j,:) = estimate_unadj(Mboot);     
  end

  se = std(beta_bs);

  % Gather coeffs
  betanace_unadj(nace,:) = beta;
  senace_unadj(nace,:) = se; 
  export_unadj = [export_unadj; export];
end

% -----------------------------------------------------
% ESTIMATION : ADJUSTED LABOR - MINCER APPROACH
% -----------------------------------------------------

export_adj = [];
for nace = naceidx 
  disp('NACE'); disp(nace);
  firmidx = M(:,9)==nace;
  Mnace = M(firmidx,:);
  Nnace = N(firmidx,:);

  [beta export] = estimate_mincer([Mnace Nnace]);

  % Boostrap standard errors
  beta_bs = [];
  for j=1:J
   disp('Bootstrap iteration'); disp(j);
    Mboot = bootstrp2([Mnace Nnace]);                            
    beta_bs(j,:) = estimate_mincer(Mboot);     
  end

  se = std(beta_bs);

  % Gather coeffs
  betanace_adj(nace,:) = beta;
  senace_adj(nace,:) = se; 
  export_adj = [export_adj; export];
end

% -----------------------------------------------------
% ESTIMATION : ADJUSTED LABOR - AVERAGE WAGE APPROACH
% -----------------------------------------------------

export_avgwage = [];
for nace = naceidx 
  disp('NACE'); disp(nace);
  firmidx = M(:,9)==nace;
  Mnace = M(firmidx,:);
  Nnace = N(firmidx,:);

  [beta export] = estimate_avgwage([Mnace Nnace]);

  % Boostrap standard errors
  beta_bs = [];
  for j=1:J
    disp('Bootstrap iteration'); disp(j);
    Mboot = bootstrp2([Mnace Nnace]);                            
    beta_bs(j,:) = estimate_avgwage(Mboot);     
  end

  se = std(beta_bs);

  % Gather coeffs
  betanace_avgwage(nace,:) = beta;
  senace_avgwage(nace,:) = se; 
  export_avgwage = [export_avgwage; export];
  
end

% ----------------------------------------------------
% ESTIMATION : ADJUSTED LABOR - GRILICHES APPROACH
% ----------------------------------------------------

[beta_gril export_gril] = estimate_griliches([M N]);

% Boostrap standard errors
beta_bs = [];
for j=1:J
  disp('Bootstrap iteration'); disp(j);    
  Mboot = bootstrp2([M N]);                            
  beta_bs(j,:) = estimate_griliches(Mboot);     
end

se_gril = std(beta_bs);

% Gather coeffs
betanace_gril(1,:) = beta_gril;
senace_gril(1,:) = se_gril;
 


%%
% -----------------------------------------
% Export data & CALCULATE TFP PREMIUM
% -----------------------------------------

betanace_unadj(betanace_unadj==0)=NaN;
betanace_adj(betanace_adj==0)=NaN;
betanace_avgwage(betanace_avgwage==0)=NaN;
betanace_gril(betanace_gril==0)=NaN;
senace_unadj(senace_unadj==0)=NaN;
senace_adj(senace_adj==0)=NaN;
senace_avgwage(senace_avgwage==0)=NaN;
senace_gril(senace_gril==0)=NaN;

disp('Unadjusted betaL, betaK, and omega coefficients');
disp((betanace_unadj));
disp('s.e.'); disp((senace_unadj));
disp('Mincer adjusted betaL, betaK, and omega coefficients');
disp((betanace_adj));
disp('s.e.'); disp((senace_adj));
disp('Average wage adjusted betaL, betaK and omega coefficients');
disp((betanace_avgwage));
disp('s.e.'); disp((senace_avgwage));
disp('Griliches adjusted betaL, betaK, quality coefficients, and omega coefficients');
disp((betanace_gril'));
disp('s.e.'); disp((senace_gril'));

pause;

% Export tfp unadjusted, tfp mincer, tfp griliches, iq griliches, tfp average wage
export = [export_unadj export_adj(:,4) export_gril(:,4) export_gril(:,5) export_avgwage(:,4)];
%format = '%d %d %d %f %f %f %f %f\n';
%fileID = fopen(strcat(workdir,'/Data/ackerberg_est2.raw'),'w');
%fprintf(fileID,format,export');
%fclose(fileID);

% Exporter premium

clear premiumunadj premiumadj ci_unadj ci_adj;
alldata.lnL = M(:,10);
alldata.Ex = M(:,8);
alldata.nace = M(:,9);
alldata.year = M(:,2);
alldata.tfp_unadj = export(:,4);
alldata.tfp_adj = export(:,5);
alldata.tfp_gril = export(:,6);
alldata.tfp_avgwage = export(:,8);

alldata.firmid = M(:,2);
alldata.lnL = M(:,10);

alldata.sex = N(:,4);
alldata.ftenure = N(:,5:7);
alldata.edu = N(:,8:11);
alldata.exp = N(:,12:15);

i = 1;
for y=min(alldata.year):max(alldata.year)
  firmidx = alldata.firmid==y;
 
  group = grp2idx(alldata.nace(firmidx));
  Dnace = dummyvar(group);
  Dnace = Dnace(:,2:end);

  n = sum(firmidx);
  X = [ones(n,1) alldata.Ex(firmidx) alldata.lnL(firmidx) Dnace];               

  [b,bint,r,rint,stats]  = regress(alldata.tfp_unadj(firmidx),X);
  premiumunadj(i,:) = [b(2)];
  ci_unadj(i,:) = [bint(2,:)];
  
  [b,bint,r,rint,stats]  = regress(alldata.tfp_adj(firmidx),X);
  premiumadj(i,:) = [b(2)];
  ci_adj(i,:) = [bint(2,:)];
  
  [b,bint,r,rint,stats]  = regress(alldata.tfp_gril(firmidx),X);
  premiumgril(i,:) = [b(2)];
  ci_gril(i,:) = [bint(2,:)];
  
  [b,bint,r,rint,stats]  = regress(alldata.tfp_avgwage(firmidx),X);
  premiumavgwage(i,:) = [b(2)];
  ci_avgwage(i,:) = [bint(2,:)];

  i = i+1;
end

disp('Exporter premium unadjusted (and confidence interval)'); disp([premiumunadj ci_unadj]);
disp('Exporter premium Mincer adjusted (and confidence interval)'); disp([premiumadj ci_adj]);
disp('Exporter premium Griliches (and confidence interval)'); disp([premiumgril ci_adj]);
disp('Exporter premium average wage (and confidence interval)'); disp([premiumavgwage ci_avgwage]);

disp('Difference in premium: Mincer vs unadj'); disp(premiumadj-premiumunadj);
disp('Difference in premium: Griliches vs unadj'); disp(premiumgril-premiumunadj);
disp('Difference in premium: Average wage vs unadj'); disp(premiumavgwage-premiumunadj);

% Plot premia by year

y1 = min(alldata.year);
y2 = max(alldata.year);
hold on;
subplot(2,1,1); line(y1:1:y2,premiumunadj);
line(y1:1:y2,premiumadj,'LineStyle','-.');
legend('TFP premium, unadj.','TFP premium, mincer');
legend('boxoff');
axis([y1 y2 -.04 .1]);
hold off;


hold on;
subplot(2,1,2); line(y1:1:y2,premiumunadj);
line(y1:1:y2,premiumgril,'LineStyle','-.');
legend('TFP premium, unadj.','TFP premium, Griliches');
legend('boxoff');
axis([y1 y2 -.04 .1]);
hold off;

%saveas(gcf,strcat(workdir,'/graphs/tfp_premia.eps'));

